
%%% Weighting the Evidence: A Rank-Dependent Model of Outdoor Recreation, June 2024
%%% This file  estimates linear and rank-dependent preferences

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IMPORTING THE DATA:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic;
Tmain=readtable('ReadyDataFluke2022.xlsx'); 
IJ=Tmain.opt_out;
ID=Tmain.ID;
C=Tmain.cost;
I=Tmain.chosen;
Z=[Tmain.release_bsb Tmain.keep_bsb Tmain.catch_scup];
R1=Tmain.catch_fluke_max;
R2=Tmain.catch_fluke_min;
R3=Tmain.keep_fluke_max; 
R4=Tmain.keep_fluke_min;
Y=[Tmain.age Tmain.income_medium Tmain.income_high Tmain.education_college Tmain.education_graduate Tmain.avidity Tmain.male];
Q=[Tmain.release_fluke_0 Tmain.release_fluke_1 Tmain.release_fluke_2 Tmain.release_fluke_3 Tmain.release_fluke_4 Tmain.release_fluke_5 Tmain.release_fluke_6 Tmain.release_fluke_7 Tmain.release_fluke_8];
G=[Tmain.keep_fluke_0 Tmain.keep_fluke_1 Tmain.keep_fluke_2 Tmain.keep_fluke_3 Tmain.keep_fluke_4 Tmain.keep_fluke_5 Tmain.keep_fluke_6 Tmain.keep_fluke_7 Tmain.keep_fluke_8];
P=[Tmain.p_keep_fluke_0 Tmain.p_keep_fluke_1 Tmain.p_keep_fluke_2 Tmain.p_keep_fluke_3 Tmain.p_keep_fluke_4 Tmain.p_keep_fluke_5 Tmain.p_keep_fluke_6 Tmain.p_keep_fluke_7 Tmain.p_keep_fluke_8];
R=horzcat(R1,R2,R3,R4); 
K=size(P,2); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAXIMIZING THE LIKELIHOOD FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rng(4560); 
a=-10;
b=10;
X0L=a+(b-a)*rand(1,14);
parameters=X0L;

%I. LINEAR MODEL:
myObjL=@(parameters)likelihoodLinear(K,Q,G,Z,P,Y,C,I,IJ,ID,parameters); 
options=optimset('MaxFunEvals',45000,'TolX',1e-35,'MaxIter',45000, 'Display','iter');
[xL,fvalL,exitflagL,outputL]=fminsearch(myObjL,X0L,options);
[XL,fvalueL,exiflagLunc, outputLunc, gradL, hessianL]=fminunc(myObjL,xL,options);
normgradL=norm(gradL); 
H_L=-hessianL; 
IM_L=-inv(H_L);
SE_L=diag(sqrt(IM_L))'; 
StatisticL=XL./SE_L; 
conditioningL=cond(H_L); 
if conditioningL<(6.7*(10^7))
condition_L='passed';
else
condition_L='failed';
end
LL0=6478.517; 
Pseudo_R2_linear=1-(fvalueL/LL0);
LL_l=fvalueL;
aic_linear=2*length(XL)+2*LL_l;
bic_linear=length(XL)*log(size(Tmain,1))+2*LL_l;

%II. CONSTRAINED RDUT MODEL (to obtain an initial value for the estimation of the T&K and Prelec models)
X0=[XL(1),XL(2),XL(3),XL(4),XL(5),rand,XL(6),XL(7),XL(8),XL(9),XL(10),XL(11),XL(12), XL(13),XL(14),1];
lb=[-Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf 1];
ub=[Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf 1];

myObj_eu=@(parameters)likelihoodGrad_RDUT_TK(K,Q,G,Z,P,Y,C,I,IJ,ID,parameters);
options=optimset('MaxFunEvals',25000,'TolX',1e-25,'MaxIter',25000, 'Display','iter','Algorithm','sqp','GradObj','on');
[X2,fval2,exiflag2, output2, lambda, grad_eu, hessian_eu]=fmincon(myObj_eu,X0,[],[],[],[],lb,ub,[],options);
H_eu=-hessian_eu; 
IM_eu=-inv(H_eu); 
SE_eu=diag(sqrt(IM_eu))'; 
Statistic_eu=X2./SE_eu; 
normgrad_eu=norm(grad_eu); 
H_eigenvalues_eu=eig(H_eu); 
conditioning_eu=cond(H_eu); 
if conditioning_eu<(6.7*(10^7))
condition_eu='passed';
else
condition_eu='failed';
end


%III. RANK-DEPENDENT UTILITY MODEL - CARA UTILITY & TVERSKY AND KAHNEMAN W(P)
X0=[X2(1),X2(2),X2(3),X2(4),X2(5),X2(6),X2(7),X2(8),X2(9),X2(10),X2(11),X2(12),X2(13),X2(14), X2(15),rand];

myObj_rdu_tk=@(parameters)likelihoodGrad_RDUT_TK(K,Q,G,Z,P,Y,C,I,IJ,ID,parameters);
options=optimset('MaxFunEvals',25000,'TolX',1e-25,'MaxIter',25000, 'Display','iter','Algorithm','quasi-newton','GradObj','on');
[X4,fval4,exiflag4, output4, grad_rdu_tk, hessian_rdu_tk]=fminunc(myObj_rdu_tk,X0,options);
H_rdu_tk=-hessian_rdu_tk; 
IM_rdu_tk=-inv(H_rdu_tk); 
SE_rdu_tk=diag(sqrt(IM_rdu_tk))'; 
Statistic_rdu_tk=X4./SE_rdu_tk; 
Statistic_rdu_tk(length(X4))=(1-X4(length(X4)))/SE_rdu_tk(length(X4)); 
normgrad_rdu_tk=norm(grad_rdu_tk); 
H_eigenvalues_rdu_tk=eig(H_rdu_tk); 
conditioning_rdu_tk=cond(H_rdu_tk); 
if conditioning_rdu_tk<(6.7*(10^7))
condition_rdu_tk='passed';
else
condition_rdu_tk='failed';
end

%Computing McFadden Pseudo-R2:
Pseudo_R2_rdu_tk=1-(fval4/LL0);
LL_rdu_tk=fval4;
aic_rdu_tk=2*length(X4)+2*LL_rdu_tk;
bic_rdu_tk=length(X4)*log(size(Tmain,1))+2*LL_rdu_tk;


%IV. RANK-DEPENDENT UTILITY MODEL - CARA UTILITY & PRELEC W(P)
X0=[X2(1),X2(2),X2(3),X2(4),X2(5),X2(6),X2(7),X2(8),X2(9),X2(10),X2(11),X2(12),X2(13),X2(14), X2(15),rand];

myObj_rdu_prelec=@(parameters)likelihoodGrad_RDUT_Prelec(K,Q,G,Z,P,Y,C,I,IJ,ID,parameters);
options=optimset('MaxFunEvals',25000,'TolX',1e-25,'MaxIter',25000, 'Display','iter','Algorithm','quasi-newton','GradObj','on');
[X6,fval6,exiflag6, output6, grad_rdu_prelec, hessian_rdu_prelec]=fminunc(myObj_rdu_prelec,X0,options);
H_rdu_prelec=-hessian_rdu_prelec; 
IM_rdu_prelec=-inv(H_rdu_prelec); 
SE_rdu_prelec=diag(sqrt(IM_rdu_prelec))'; 
Statistic_rdu_prelec=X6./SE_rdu_prelec; 
Statistic_rdu_prelec(length(X6))=(1-X6(length(X6)))/SE_rdu_prelec(length(X6)); %t-test again the null that, in the  weighting function, gamma=1.
normgrad_rdu_prelec=norm(grad_rdu_prelec); 
H_eigenvalues_rdu_prelec=eig(H_rdu_prelec);
conditioning_rdu_prelec=cond(H_rdu_prelec); 
if conditioning_rdu_prelec<(6.7*(10^7))
condition_rdu_prelec='passed';
else
condition_rdu_prelec='failed';
end

%Computing McFadden Pseudo-R2:
Pseudo_R2_rdu_prelec=1-(fval6/LL0);
LL_rdu_prelec=fval6;
aic_rdu_prelec=2*length(X6)+2*LL_rdu_prelec;
bic_rdu_prelec=length(X6)*log(size(Tmain,1))+2*LL_rdu_prelec;


%%%%%%%%%%%%%%%%SAVING THE OUTPUTS OF BASELINE MODELS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dataL=horzcat(vertcat("LINEAR","-","-"),vertcat("Parameter","SE","test stat"),vertcat(XL,SE_L,StatisticL),vertcat("LogL","AIC","BIC"),vertcat(-LL_l,aic_linear,bic_linear));
data_rdu_tk=horzcat(vertcat("T&K","-","-"),vertcat("Parameter","SE","test stat"),vertcat(X4,nonzeros(SE_rdu_tk)',Statistic_rdu_tk),vertcat("LogL","AIC","BIC"),vertcat(-LL_rdu_tk,aic_rdu_tk,bic_rdu_tk));
data_rdu_prelec=horzcat(vertcat("PRELEC","-","-"),vertcat("Parameter","SE","test stat"),vertcat(X6,nonzeros(SE_rdu_prelec)',Statistic_rdu_prelec),vertcat("LogL","AIC","BIC"),vertcat(-LL_rdu_prelec,aic_rdu_prelec,bic_rdu_prelec));
V1=array2table([max(ID),length(unique(Tmain.qtid));NaN,NaN;NaN,NaN],"VariableNames",["N. choice occasions","N. anglers"]);
V2=array2table([NaN,NaN;NaN,NaN;NaN,NaN],"VariableNames",["N. choice occasions","N. anglers"]);
Table_Results_linear=array2table(dataL,"VariableNames",["Model","Variable","SFrelease", "SFkeep", "BSBrelease", "BSBkeep","Scupcatch","Cost", "Constant Opt_out" ,"Age", "Income_medium", "Income_high", "Education_college","Education_graduate",'Avidity','Male','Fit',"N. choice occasions"]);
Table_Results_rdu_tk=array2table(data_rdu_tk,"VariableNames",["Model","Variable","SFrelease", "SFkeep", "BSBrelease", "BSBkeep","Scupcatch","Risk Aversion", "Cost", "Constant Opt_out" ,"Age", "Income_medium", "Income_high", "Education_college","Education_graduate","Avidity","Male","Probability weight","Fit","-","N. choice occasions"]);
Table_Results_rdu_tk=[Table_Results_rdu_tk V1];
Table_Results_rdu_prelec=array2table(data_rdu_prelec,"VariableNames",["Model","Variable","SFrelease", "SFkeep", "BSBrelease", "BSBkeep","Scupcatch","Risk Aversion", "Cost", "Constant Opt_out" ,"Age", "Income_medium", "Income_high", "Education_college","Education_graduate","Avidity","Male","Probability weight","Fit","-","N. choice occasions"]);
Table_Results_rdu_prelec=[Table_Results_rdu_prelec V2];
writetable(vertcat(Table_Results_rdu_tk,Table_Results_rdu_prelec),'N:\Faculty\JHolzer\ado\McConnell_Holzer_RDUT\Matlab_code\estimation\Table2.csv')

save('paramlinear.txt','XL','-ascii');  save('parameters_rdu_tk.txt','X4','-ascii');  save('parameters_rdu_prelec.txt','X6','-ascii'); 
save('sigmaL.txt', 'IM_L', '-ascii'); save('sigma_rdu_tk.txt', 'IM_rdu_tk', '-ascii'); save('sigma_rdu_prelec.txt', 'IM_rdu_prelec', '-ascii');   
save('ModelEstimates'); 

toc;

toc;